####################################### SINGLE EXPOSURE ##################################################


#### Opening the data ####

# You need to put your address
aids<-read.table("C:/Users/laupaz/Desktop/Teaching & visiting/Causal inference (Milano Bicocca) Dicembre 2018/Lab/single.txt",header=TRUE)
# Visualizing the first 10 rows
aids[1:10,]


####### Unadjasted analysis in R of the relation among A and Y ########

chisq.test(x=aids$A,y=aids$Y)

####### Adjasted analysis in R of the relation among A and Y########

# Fitting the model in R without interaction:
summary(glm(formula=Y~A+L,family=binomial,data=aids))

# Adding an interaction term between A and L:
summary(glm(formula=Y~A+L+A*L,family=binomial,data=aids))


###### Marginal effect through outcome model: ML ######

# STEP 1 
# Fit a regression model for the outcome
fit0<-glm(formula=Y~A+L+A*L,family=binomial,data=aids)

# STEP 2 
# Replace the factual exposure level with a, for each individual 
# (for a=0)
aids0=aids
aids0$A=0

# STEP 3 
# Estimate Pr(Y = 1|A = a; L) for each individual (for a=0)
# Installing package 
library("stats", lib.loc="~/R/R-3.3.0/library")
pred0<-predict(object=fit0,newdata=aids0, type="respons")

# STEP 4 
# Average these estimates over all individuals to obtain an estimate 
# of Pr(Ya = 1) (for a=0)
p0<-mean(pred0)
p0

# STEP 1 
# Fit a regression model for the outcome
fit1<-glm(formula=Y~A+L+A*L,family=binomial,data=aids)

# STEP 2 
# Replace the factual exposure level with a, for each individual 
# (for a=1)
aids1=aids
aids1$A=1

# STEP 3 
# Estimate Pr(Y = 1|A = a; L) for each individual (for a=1)
pred1<-predict(object=fit1,newdata=aids1,type="respons")

# STEP 4 
# Average these estimates over all individuals to obtain an estimate 
# of Pr(Ya = 1) (for a=1)
p1<-mean(pred1)
p1

# Marginal causal log odds ratio
Marg_clor<-log((p1/(1-p1))/(p1/(1-p0)))

####### Marginal effect through exposure model: IPW #######

# STEP 1 Fit a regression model for the exposure
fit<-glm(formula=A~L,family=binomial,data=aids)

# STEP 2 For each subject, use the fitted exposure model to
#estimate a subject-specific weight
pred<-predict(object=fit,type="respons")
w<-1/(aids$A*pred+(1-aids$A)*(1-pred))

# STEP 3 Use Pr(Y = 1|A = a) in the weighted sample as an
# estimate of Pr(Ya = 1)for a=0
p0<-weighted.mean(x=aids$Y[aids$A==0],w=w[aids$A==0])
p0

# STEP 3 Use Pr(Y = 1|A = a) in the weighted sample as an
# estimate of Pr(Ya = 1
p1<-weighted.mean(x=aids$Y[aids$A==1],w=w[aids$A==1])
p1

# Marginal causal log odds ratio
Marg_clor<-log((p1/(1-p1))/(p1/(1-p0)))


# Repeat the estimation of marginal effects adding a square term for CD$ (L)
# You can repeat all the steps above adding I(L^2) to both the outcome and exposure models